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SUMMARY 


In the present study, a scheme capable of solving very fast and robust 
complex nonlinear systems of equations is presented. The Block Adaptive 
Multigrid (BAM) solution method offers multigrid acceleration and adaptive 
grid refinement based on the prediction of the solution error. The proposed 
solution method was used with an implicit upwind Euler solver for the solution 
of complex transonic flows around airfoils. Very fast results were obtained 
(18-fold acceleration of the solution) using one fourth of the volumes of a 
global grid with the same solution accuracy for two test cases. 


INTRODUCTION 


Although multigrid methods were introduced as grid adaptation techniques 
they have been established only as fast and efficient solvers for large scale 
computational problems. Up until today only a few adaptive multigrid schemes have 
been presented, e.g. the Multilevel Adaptive Technique MLAT (ref. 1), the Fast 
Adaptive Composite grid method FAC (ref. 2) and others (ref. 3), but the domain 
of applications has been mainly restricted to the solution of elliptic type 

equations. Regarding the development of adaptive schemes for hyperbolic 
systems of equations, only a few attempts have been made to take advantage of the 
favourable multigrid concept for the acceleration of the solution. On the 

other hand great advantages have been pointed out for the use of the 

truncation error prediction as a reliable error sensor for grid adaptation, 
though few studies exhibited numerical proofs (ref. 4,5). 

The present study, a dynamically grid adaptive method, namely the Block 
Adaptive Multigrid (BAM) method, is presented, incorporating a reliable device 

for the prediction of the error and a composite multigrid solver. The method 
is based on a Full Multigrid scheme: starting from an acceptable coarse mesh 

the solution creates finer grid patches with block grid refinement. The 

refined regions are grouped into rectangular blocks defining a composite 

structure which is totally handled by the multigrid method. In this way an 

adapted non uniform domain is decomposed into regular subdomains solved with 
common solvers. Further, the communication between blocks and the 

1 This work was partially supported by CEC/ Brite-Euram contract AERO-0018C. 
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parallelization of the BAM method are based on domain decomposition ideas. For 
the integration and relaxation of the time marching Euler equations, an 
unfactored implicit upwind finite volume scheme is employed. The proposed 
method is tested for two complicated transonic inviscid cases for two 
different airfoils. Using the proposed method stable and accurate results are 
obtained in a small number of work units (18-fold acceleration with 4.5 times 
fewer volumes). 


FINITE VOLUME DISCRETIZATION 


The general inviscid flow is described by the Euler equations that can be 
solved using the very popular time-marching conservative formulation. For the 
two dimensional case, conservation laws are used with body- fitted coordinates 
& fl: 
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The steady state solution is found when the time derivative of the solution 
vector becomes negligible. The solution vector and the fluxes normal to 
£= const, and q = const, faces are given respectively: 


U - J • U and E - J • (E ^+F £ y ) , F = J • (E + Fr^) (2.2) 

At the Cartesian coordinate system the corresponding solution vector and the 
inviscid fluxes are given by: 
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In equation (2.3) e is the total energy (e = (y^) +^-p[u 2 +v 2 ]), p and p are 

the pressure and the density respectively and J is the Jacobian of the inverse 
mapping. 


For the discretisation of equation (2.1) a cell-centered finite volume 
method is used. For the pseudo-time evolution a Newton linearization scheme is 
adopted which, being an implicit scheme, allows high CFL numbers (100-200) to 
be used with local time stepping (different At for each volume): 

-ST- * < A " AU \ + < B " AU ) n - - (E"+ 9 

Or else: 

L n - AU = ^At* 1 + Aj) + B”j • AU = - Res n (U n ) 

Where AU is the correction of the solution vector, A and B are the Jacobians 


(2.4) 

(2.5) 
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( 2 . 6 ) 


of the fluxes E and F respectively for the time level n: 
u" + 1 - u" ♦ AU . A" = (-§£-)" and B" = (-§£-)" 

Upwind differencing of the flux vectors is used to achieve a diagonal 
dominant system of equations. Thus, using a symmetric collective point Gauss- 

Seidel relaxation scheme very good smoothing properties are attained for the 

multigrid calculations. For the flux calculations a linear locally one- 

dimensional Riemann solver (Godunov- type) is employed thus, the homogeneous 
property of the Euler fluxes [7] is guaranteed. The mean values of the 

conservative variables at both sides of the faces are used as input variables 
for the Riemann solver. For the calculation of the fluxes E and F the 

conservative variables are extrapolated using an upwind characteristic 

variable interpolation method (MUSCL- type). The interpolation scheme uses two 
volumes from both sides of each face. The accuracy of the scheme raises, up to 
third order depending on the sign of the local eigenvalues of the Jacobians A 
and B. The local accuracy of the finite volume method is sensor- controlled so 
the monotonic behaviour of the solution is guaranteed. Boundary conditions are 
required for both sides of eq.(2.4), so for the RHS of eq.(2.4) 
characteristic boundary conditions are extracted from the Riema.nn solutions at 
the wall and at free surfaces. For the LHS of eq.(2.4) simple boundary 
conditions are prescribed on the AU in phantom shells. 


THE BLOCK ADAPTIVE MULTIGRID METHOD 


The Block Adaptive Multigrid (BAM) method is composed of three main 
parts: the fast nonlinear multigrid solver, the truncation error approximation 
for the prediction of the solution error and the block composite grid solver. 
Additionally an efficient solution strategy is required in order to achieve 
fast and robust accurate solutions. 


Multigrid Implementation 


Concerning the multigrid method the Full Approximation Scheme (FAS) is 
employed as it is better than the Correction Scheme (CS) for Newton 
linearisation schemes. FAS has the major advantage in that it operates with the 
same solution vector of the initial algorithm so it is best suited for the 
solution of composite grid structures. The efficiency and the performance of 
the multigrid implementation are maintained adopting the "alternate point of 
view" of the FAS (ref. 1). Following this approach the finer grid levels are 
considered as devices to increase the spatial accuracy of the solution whereas 
the coarser grid levels are devices to accelerate the solution. The 
formulation of the solution is independent of the grid level (coarse or fine) 
and the type of grid (local or global) by simply adding^ to the RHS of eq.(2.5) 
the appropriate fine-to-coarse defect correction (x). For the enumeration of 
the multigrid levels the classical mode is adopted. Thus, the current grid 
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level is denoted by n, its next finer level by n + 1 and its corresponding 
finest grid level by m (m s n s 1) (1 is the coarsest grid level of the 
domain). The solution formulation for the current grid level n is given as: 


n+l 


L • AU = - Res (U ) + x 1 

n n n v n' n 

considering that the fine-to-coarse defect correction is: 

x "n 1= *l'\ L „ + r AU n + i 1 - L • (I n+1 AU .) and 
n n [ n+l n+l J n v n n+l y 

Because one multigrid cycle equals one time step the time scale is not 
considered. 
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Because of the applied cell centered finite volume scheme the cellwise 
coarsening is adopted; each coarse volume is constructed by four consequently 
finer ones. The cellwise coarsening maintains the outer edges of the volumes 
(straight implementation of conservation laws) but the coarse grid centers are 
not a subset of the fine ones. Therefore, two different restriction operators 
are required. The restriction operator (I) for the physical variables is the 
simple average over the four fine volumes: 
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In contrast, the restriction operator (E) for the generalized residuals, Res 
and x, is the summation of the residuals of the corresponding fine volumes. The 
fluxes of the inner common fine grid faces are canceled, thus flux 
conservation at the coarser grid levels is maintained: 


Res = E n + 1 Res 

n n n+l 



(3.4) 


For the reverse direction of the multigrid cycle (coarse-to- fine direction) 
neither Euler solutions nor relaxation sweeps are required. Therefore, AU 
variables are stored for all grid levels and only these are prolongated from 
the coarse-to-fine direction using the standard FAS prolongation formulation: 
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For the prolongation operator (II) simple injection is adopted, i.e.: 


(3.5) 
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Due to the composite grid structure, complicated interpolation schemes would 
have increased considerably the programming complexity with minor advantages. 
For the present multigrid implementation the V-cycle is applied with the 
improvement of increasing the number of relaxation sweeps as coarser levels 
are processed. 


Solution Error Approximation 


To determine the erroneous regions of the computational domain where 



increase of the accuracy is required, a reliable error indicator . 1S req “ ire ‘ 
Two approaches essentially exist: the first one is the physically based 

information of the problem, i.e. solution gradients while the second one 
numerically motivated, is the evaluation of the discretization error. Th 

former approach may implicate the refinement process to reduce errors that 
have no influence on the global solution. In contrast, the evaluation of 
truncation error approach indicates errors which can be confronted by 
refinement procedure. On the basis of the Richardson extrapolation 
estimation of the truncation error is formulated as: 

T = Q(h)*u ( 3 - 6 ) 

n 

Whereas for a uniform Cartesian mesh holds: 


T - c»(h) p 


(3.7) 


In eauations (3 6) (3.7) Q is the differential operator, u is the physically 

^orrS So n is the mesh size of the finest grid level, p is the local 
order of accuracy and c is an unknown grid independent to»jOuttW ' 
nhvsical interpretation of the truncation error and the Richardson 

extrapolation concept, the difference of residuals between 

levels is used for the truncation error calculation. This error sensor 
provides a reliable local estimation of the solution error. ,AM° u gh the 
fine-to-coarse defect correction of eq.(3.2) is directly provided J 

multigrid solution it was proved to be an unreliable error indicator. 

Fortunately the use of the Res(U) operator instead of the L°AU operator gves 
very good results. The explanation is that due to the Newton lineariza 
scheme the Res(U) differential operator is insensitive to relaxation errors 
maintaining the accuracy of the solution. Therefore, the solution error 

evaluation ® for the grid level m-1 (m is the current local finest grid level), 

is given by: 

(3.8) 


T = Z m Res (U ) - Res TI .U ) 

‘m-l — m - 1 "m m- 1 m-1 m m m-1 m-1 m 

For a totally converged solution equation (3.8) reduces to. 
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As this truncation error sensor is a vector, a reduction norm to a single 
value is required. At the present study the Euclidean norm has been adopted 
because it shows similar distribution to the pressure error of the solution. 
The proposed error sensor requires additional work of only one fourth of 
simple r flux calculation and it does not demand totally converged solution 
(Res (U )*0) although it converges from the early time steps to the stea y 


state. 


The prediction of the solution error for the new finer grid level m-1 can 
be computed only if the assumption of a uniform Cartesian mesh is adop 

(h = h 12). Thus following equation (3.7) we get: 
v m+1 m 
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For the determination of the order of accuracy (p) of eq.(3.10) which 
varies from one to two depending on the flow features, the sensing functions 
that are used to control monotonicity at the integration routine are also use 
to calculate the local accuracy of the scheme. Unfortunately, equation (3.10) 
holds only for Cartesian grids, thus using this equation in arbitrary gri s 
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errors are expected to the predicted error levels. 


Composite Grid Structure and Solution Procedure 


. . For the optimal approximation of the most accurate solution within the 
minimum amount of work, several grid adaptive techniques and structures have 
been developed. The composite grid structure has many advantages by enabling 
the solution of a globally non uniform grid using a union of locally uniform 
subgrids (Mocks) (ref. 2). Subgrid uniformity is essential to assure 

multigrid efficiency using simple solution routines (similar to the single 
grid solver). Moreover, significant simplifications to the data structure and 
to the interface manipulation are attained when the blocks are restricted to 
rectangles in the computational domain and the grid refinement ratio is 2 for 
each refined direction (ref. 5). 


. a composite multigrid method the major problem is the requirement of 
special manipulation at the artificial boundaries to maintain the accuracy of 
the global solution. To suppress errors inconsistent with the solution method 

at the artificial boundaries (two fine volumes share the same edge with a 

coarse one) certain special boundary conditions must hold. Namely: accuracy 
preservation of the integration scheme, flux conservation and flux- splitting 
compatibility with respect to the global grid solver. To preserve the accuracy 
of the solution across an artificial boundary, the fluxes of the domain should 
be calculated at the block’s finest grid level. For these calculations the 
standard integration routine has been employed. The basic idea is to properly 
construct false fine volumes at the coarse grid boundaries using all the 

available information near the interfaces. The modified scheme at the 
intergrid boundaries is depicted in fig.l and fig.2 for the one and two 

dimensional analogs, respectively. Following these analogs the flux 
calculations should be started from Block 1 including the boundary interfaces 
For die evaluation of the virtual volumes F3 and F4 (fig.l) the coarse volumes 
Cl, C2 and the fine ones FI, F2 are considered adopting the same MUSCL- type 
extrapolation scheme used by the integration routine. At the two dimensional 
analog, the same couple of coarse volumes is used with the corresponding fine 
volumes (fig.2). In this way, the order of accuracy of the global solution 

scheme is maintained. After calculating the fluxes of the finest subgrid 

1 in the boundary fluxes of the neighboring coarser subgrids 

(Block 2) can be calculated explicitly using conservation of fluxes. According 
to the multigrid restriction operator for the residuals, flux conservation 
across the artificial boundaries is achieved by addressing the summation of 
the fluxes of the fine volume faces to their adjacent coarse volume face (fig. 

With respect to the relaxation procedure, the modifications to the flux 
vector splitting and the relaxation scheme at the subgrid interfaces are 

totally handled by the multigrid algorithm. Adopting the "horizontal" 

communication mode among subgrids, the fine subgrids (i.e. block 1) are 

relaxed at the first grid level (fig.l). At the next multigrid level the 

coarse subgrids (i.e. block 2) are relaxed together with the restricted fine 

ones, while the block structure of the composite grid is preserved throughout 

the multigrid cycle. 
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Due to the block structure, the composite grid has the advantage of a 
straightforward implementation of vectorization and paraUelization. This can 
be achieved when subgrids are considered totally independent from each other 
introducing ideas from the domain decomposition theory. Concerning the 
"vertical" communication mode (ref. 11), each subgrid is solved and relaxed 
independently for the complete multigrid cycle whereas data exchange amog 
subgrids is permitted only at the start and/or the end of each cycle. TTie 

computational domain decomposition (different from the composite bl 
structure) (ref. 8) should be performed according to load balancing and 
vectorization criteria taking into consideration the hardware configuration. 

The "vertical" communication mode is preferred from the "horizontal mo e w en 
we deal with parallel processing and is used even though the latter mode has 

better convergence rates. Using the "vertical" mode the idle and communication 
time among processors can be considerably decreased. 

With respect to the dynamicaUy adaptive multigrid strategy a modified 

Full Multigrid scheme is applied, starting with a global coarse grid o 
acceptable grid resolution. After convergence (or after a fixed amount of 
work) at the current grid, the truncation error is calculated and the sohition 

error is predicted. In regions where the prediction of the error is above a 
threshold the corresponding volumes are flagged and grouped into rectangular 
blocks. Afterwards, the domain is decomposed to the appropriate blocks where 
only those which contain the flagged volumes are refined to the next grid 
level injecting the coarse grid solution to the refined grid. TTie refinement 
procedure continues until the entire computational domain has local truncation 
errors below a given threshold. Clearly this strategy has the benefit of a 
continuous iterative procedure without wasting CPU-time on calculations that 
will not be used after the mesh refinement. Taking advantage of the most 
accurate available solution the proposed scheme converges straight to the most 
efficient solution. 


RESULTS AND DISCUSSION 


In order to verify the accuracy and to validate the efficiency of the 

proposed method two transonic inviscid test cases were investigated. The first 
case is a NACA- 0012 airfoil for Mach 0.80 at angle 1.25 degrees while the 

second case is a RAE- 2822 airfoil for 0.73 Mach at 2.79 degrees. A Work Unit 
fWU) is defined by the CPU-time required for a global finest grid relaxation 

sweep ^ in lexicographic order while for a single grid running one time step 
costs four work units. 

For both test cases the grid adaptation procedure as described above was 
applied. Starting point is two global multigrid levels with 64x14 volumes at 
the finest grid level. The user supplies only the maximum number of the 
additional grid levels which for both test cases were the same: two refinement 
grid levels The truncation error threshold was defined explicitly targeting 

fo a three- fold reduction of the initial error levels. For the convergence 
criterion the Euclidean norm of the correction vector was employed. 

For the first test case (NACA-0012) a comparison between the computed and 
the predicted truncation error contours is given in figures 3 and 4 - re *P" 

Some differences, not crucial, are due to the incorrect evaluation of the 
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local accuracy of the solution (p). For the same grid level the actual 
pressure error contours and the entropy contours in figures 5 and 6 are 
included. A comparison between the truncation errors and the pressure errors 
shows their similar distributions. Following the proposed method very accurate 
results were obtained in comparison to the global solutions as depicted in 
fig.7 In fig. 7 the Mach- distributions along the airfoil are depicted for two 

global grids (256x56 and 128x28) and one composite grid sharing both grid 
levels. A comparison of these Mach distributions shows that at regions 

* , the same mesh size the solutions between a global and a local refinement 

coincide. Additionally, in fig.8 the Mach contours together with the composite 
grid^ are given. In fig.9 the efficiency of the BAM method with respect to the 
single grid and the global multigrid schemes is clearly indicated. A 19-fold 

and 4-fold acceleration were achieved with respect to the single grid and to 
the multigrid cases, respectively. The final grid adaptive solution requires 

4.5 times fewer volumes (from 14336 to 3200) for practically the same accuracy 
with a globally refined grid (0.35% discrepancy of the computed Cl). 

Similar efficiency was achieved for the second test case. The domain 

decomposition into subgrids took place after 50 time steps spent on coarser 

grid levels. Using the same truncation error threshold as in the previous test 

case, a 17-fold acceleration was achieved with respect to the single grid and 

a 4.7- fold reduction of the volumes (from 14336 to 3056) for practically the 
same accuracy (Cl discrepancy is 0.2 %). The convergence histories of the 
error reduction and the lift coefficient are shown in fig.10 while in fig.ll 
the final composite grid together with the isomach contours are depicted It 
is important that throughout the solution process the multigrid convergence 
rates were maintained while the overhead for the interface computations was 

negligible, i.e. only 2 % for a nine block structure with respect to an 

equivalent global grid. 

As no parallel machine was currently available a simulation of the 

parallelization for the procedure has been attempted in order to evaluate the 
performance of the two different communication modes for the BAM method. To 
achieve this, the data exchange was properly adjusted among subgrids according 
to each communication mode. For the "horizontal" mode (semi-parallel) and the 

"vertical" mode (parallel) the convergence histories are shown in fig. 12. A 
comparison between the parallel modes and the sequential BAM method shows only 
a small reduction of the convergence rates for the parallel ones. What is 
further required for a parallel implementation is a computational domain 
decomposition of the composite structure based on specific load balancing 
criteria in order to reduce the idle time among processors. & 


CONCLUSIONS 


The great advantages of the Block Adaptive Multigrid (BAM) method were 
exhibited. The incorporation of numerous efficient schemes into the BAM method 
makes the ultimate target of solving complex problems in just a few work units 
feasible. At the same time robustness, simplicity and accuracy of the single 
grid code were maintained at the new method. The extension to viscous and 
hypersonic three dimensional problems is straight forward though semi 
coarsening multigrid can be also included. On the other hand, in order to 
improve the adaptation capabilities, a moving grid point scheme should also be 


considered as the grid alignment towards certain flow features is essential in 
some problems in combination with the present grid refinement procedure. 

Although the basic features of the BAM method have been determined and 

the local order of accuracy of the solution in such a way that the prediction 
f thp solution error would be more accurate. Additionally, the implementation 
of the BAM method to other solution algorithms and equations is ° res ®®|' ‘ 

°he BAM method was designed in the general concept of the finite volume 

method. 
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Figure 1. One dimensional analog of the multigrid coarsening and the 
fictitious volumes F3.F4 for the intergrid flux calcultations. 


• 

to b 

4 


F5 

F6 

* y * 

• Cl 

4 <4 



X.C2 

M 'A 

0 

to b 

FI 

F2 

^ 



Figure 2. Two dimensional analog of the special variable interpolation procedure at the 
integrid boundaries. 
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Figure 3. Computed truncation error contours for the finer grid level, CASE 1. 



Figure 4. Predicted truncation error contours for the finer grid level, CASE 1. 
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Figure 7. Mach distribution along the airfoil for a global fine 
a global coarser and an adaptive composite grid (CAob I;. 
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Figure 8. The adaptive composite grid and the corresponding Mach contours (CASE 1) 






478 


LIFT COEFFICIENT 





Figure 11. The adaptive composite grid and the corresponding Mach contours (CASE 2). 
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Figure 12. Convergence histories of the lift coefficient and the logarithmic Euclidean error 
reduction of the B.A.M. method (CASE 2) for sequential processing, “parallel” horizontal 
and “parallel” vertical schemes. 
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